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We find a statistical mechanism that can adjust orientations of intracellular filaments to cell 
geometry in absence of organizing centers. The effect is based on random and isotropic filament 
(de-) polymerization dynamics and is independent of filament interactions and explicit regulation. It 
can be understood by an analogy to electrostatics and appears to be induced by the confining bound- 
aries; for periodic boundary conditions no orient at ional bias emerges. Including active transport 
of particles, the model reproduces experimental observations of vesicle accumulations in transected 
axons. 



In living cells, cargo like nutrients, proteins or or- 
ganelles, has to be carried to distinct locations in the 
cell. Many of these objects are transported by motor 
proteins (kinesin, dynein and myosin) that bind cargo 
and move along polarized intracellular filaments (micro- 
tubules, actin filaments) in a directed manner [1]. The 
direction of movement is determined by the polarity of 
the filament (plus- and a minus-end) and the motor pro- 
tein species. 

For particular cell functions, accumulations of cargo 
and concentration gradients are needed at distinct lo- 
cations of the cell. For that purpose the arrangement 
of filaments, which determines the distribution of cargo, 
has to be organized properly. In some cases filaments are 
organized by auxiliary components to form globally po- 
larized structures (e.g. the microtubule organizing center 
[2]). Next to these highly regulated mechanisms, mo- 
tor driven transport processes take place also on less 
organized networks. In this article we consider a sce- 
nario where the initially ordered network of microtubules 
(MTs) in a neuronal axon is destroyed after transection 
of the axon. Experimental observations show that after 
transection, vesicles agglomerate adjusting to the geom- 
etry of the cut axon [3, 4]. This vesicle agglomeration 
is important for the formation of a growth cone that is 
needed to recover neuronal connections. It was supposed 
that a re-organization of the microtubule network creates 
effective kinetic traps for the vesicles. 

Mechanisms for self-organization of intracellular fila- 
ments were proposed in previous works, relying on extrin- 
sic chemical gradients induced by external signals [5, 6], 
or spontaneous structure formation driven by the effec- 
tive interactions of filaments mediated by molecular mo- 
tors [7-11]. 

Here, we use a complementary approach: We ask the 
question how structures self-organize without relying on 
explicit regulation mechanisms (e.g. by pre-established 
chemical gradients) or interactions, and discuss if they 
are necessary to explain the alignment of filament ori- 
entations and subsequent vesicle agglomerations in finite 
volumes. We consider explicitly a simple homogeneous 
and isotropic stochastic model for filament growth dy- 



namics that neglects filament-filament interactions [12]. 
The model considers basic filament nucleation and (de-) 
polymerization dynamics as considered in previous mod- 
els (e.g. [13]), without an a priori preferred direction or 
position of filaments nor any gradients. However, the net- 
work and particle dynamics are assumed to be confined 
to a given geometry corresponding to the finite cell vol- 
ume. We show that the confining boundary conditions 
in fact induce an orientation of the filament ensemble. 
We also present a linear theory for filament network evo- 
lution which is analogous to electrostatics and able to 
explain the observed orientation of filaments. 

Our model is defined analogue to [12]. We apply 
stochastic dynamics in continuous space for filament and 
particle dynamics. The model rules capture the basics 
of filament nucleation and (de-) polymerization dynamics 
[1]: Filaments are nucleated with spatially homogeneous 
rate oj n p mon p nuc and isotropic random orientation. Here 
the values pm 0n and p nuc denote the (conserved) quantity 
of monomers and nucleation seeds in the cytosol. Fila- 
ments can grow by including segments of length d s at 
their plus-end with rate oj g p mon and shrink by dissoci- 
ating a segment at the minus- end with rate uj s . After 
nucleation, the minus-end remains capped for some time 
such that depolymerization at the minus-end is not pos- 
sible. This restriction is removed with uncapping rate uo u 
(which is non-zero only for actin filaments). 

Particles are implemented as hard-core spheres with 
radius r p exhibiting mutual steric exclusion. Detached 
particles perform a continuous random walk in space. If 
they are within the binding range of a filament d\>, they 
can attach to the filament (rate uj a ) and perform directed 
motion until they detach with rate ujd- We consider two 
species of particles: plus-particles move towards the plus- 
end of the filament with rate p per segment, while minus- 
particles move to the minus ends. The model dynamics 
are illustrated in Fig. 1. For a more detailed definition 
we refer to the supplementary material and [12]. 

We want to apply the model to the experimental setup 
studied in [3]: Healthy axons usually contain bundles of 
axially orientated microtubules with plus ends pointing 
towards the synapse. In the experiments, axons were 
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FIG. 1. Illustration of model dynamics, (a) Filaments are 
implemented as sequences of linear segments. These are polar- 
ized, with a plus-end where segments are created to elongate, 
and a minus-end where segments dissociate causing shrink- 
ing, (b) Dark gray discs are detached particles, black discs 
represent particles attached to filaments. 
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FIG. 2. Comparison between experimental results of Erez 
et al. [3] (a) and model results (b). The same color-coding 
is used: anterograde (plus-) particles are green, retrograde 
(minus-) particles are red. The shape of the model axon is 
based on the experiment (bar=15 fim). The deviations from 
the cylindrical shape have only little influence on the par- 
ticle distribution, hence the accumulation is not result of a 
(diffusive) bottleneck (see supplementary material). 



transected leading to dissolution of microtubules near the 
location of the cut, due to calcium influx (right part of 
Fig. 2(a)). In that region, microtubules subsequently 
reassemble randomly in absence of an organizing center. 
Non-destroyed MTs remain and can grow again. One 
observes an accumulation of plus- (anterograde) vesicles 
at the tip and at some distance from the tip as well as 
minus- (retrograde) vesicles between these locations. 

In the model we implement this scenario by introduc- 
ing a confining volume that is based on the cylindrical 
shape of an axon. We use an idealized cylindrical shape 
as well as a shape similar to the contour of the axon used 
in the experiment of Erez et al. [3]. The finite volume 
leads to a straightforward modification of the filament 
dynamics: Filaments are not allowed to grow or nucleate 
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FIG. 3. Comparison between experimental [3] (top) and sim- 
ulated particle distribution (bottom), z = cylinder axis coor- 
dinate. 



outside the volume (see Fig. 2(b)). Starting from paral- 
lel bundles of microtubules, they get dissolved (transec- 
tion) within some distance from the tip (right half of the 
figure) and new filaments reassemble randomly without 
preferred direction, according to the rules defined above. 
Moreover, we assume that filament growth is inhibited by 
the steric hard-core particles such that filaments are not 
allowed to grow through a particle. In addition, we imple- 
ment microtubule catastrophe-rescue dynamics following 
the model in [13]. At the left boundary, plus particles 
enter the system with rate ql\ where also minus particles 
exit the system with rate /?2- At the other boundaries 
(representing the membrane), plus particles exit (exocy- 
tosis) with rate j3\ and minus particles enter the system 
(endocytosis) with rate a^- 

In Fig. 2(b) a particle configuration 24000 time steps 
after transection (corresponding to about 4 min in real 
time) is displayed. We used the biologically motivated 
parameters as discussed in the supplementary material. 
Starting with parallel filaments, they are dissolved in the 
right part of the system after 7000 time steps, whereupon 
random nucleation and growth of filaments, according to 
the dynamics discussed above, takes place in this region. 

The comparison between experimental and model re- 
sults shows that the model is able to reproduce the struc- 
ture of the plus particle traps as well in the middle of the 
system and at the tip, but it underestimates its strength 
(Fig. 3(b)). 

In order to understand the emergence of particle ac- 
cumulations, we develop a linear theory for the network 
structure. The orientations of filament segments will be 
denoted by a unit vector d describing the direction of the 
plus-end. We quantify the structure of the network by a 
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vector field 



F(x) := ( lim 



V ^ 

fil. in v 



d(x) =p s (d(x)) (1) 



which we call filament field. Here, v is a volume con- 
taining the point x. The sum includes all filament seg- 
ments in v, p s denotes the density of filament segments 
at the point x and (• • • ) denotes the ensemble average. 
The direction of F represents the average orientation of 
filaments, whereas its amplitude is determined by the 
density of filaments. F(x) results from a linear superpo- 
sition of filament orientations at the point x which allows 
to express it by a Green's function as 



F(x)= / G(x-x / )p n (x / )^ /3 
Jv 



(2) 



where G(x — x') is the field if nucleation of filaments 
would only be possible in a single point x', and p n (x) 
denotes the relative nucleation rate. We consider the 
dynamics to be spatially homogeneous within the volume 

V where filaments can nucleate, hence p n ( x ) = 1 within 

V and p n (x) = outside [14]. However, (2) can also be 
applied, when nucleation dynamics vary spatially. 

Since the orientation of filaments always points away 
from the nucleation point, the Green's function G(x), 
which considers only nucleation in the origin, must have 
a radial structure G(x) = G(x)e r . In fact, G(x) is 
the probability that a filament is within a volume unit 
dx 3 around x. Since filaments are not correlated, this 
can be expressed by the probability P(x) that the fil- 
ament is directed towards x and, assuming ergodicity, 
the average time r(x) the filament is present in dx 3 : 
G(x) = P(x)r(r). Because a filament grows isotropically 
in arbitrary direction, the probability that it passes the 
volume dx 3 located at distance r = |x| from the origin is 
given by P(x) = P(r) oc 1/r 2 . In the stationary state, 
the net polymerization rate 0J g p mon must equal the de- 
polymerization rate uj s due to conservation of monomers. 
If we neglect shrinking at the plus end, a filament can 
only disappear at a given point by minus-end depolymer- 
ization. The average filament length I is independent of 
the distance r and hence the time the filament is present 
at x is r = l/uo s which is also independent of r. 

After all, the overall value of the filament field is 



G(r) oc 



(3) 



This form corresponds to the Green's function of elec- 
trostatics (Coulomb law). As a result, the filament field 
inside a volume V with homogeneous nucleation rate cor- 
responds to a electrostatic field of a homogeneous charge 
distribution within V. 

In order to test these results quantitatively, we con- 
sider more simple boundary conditions than in the pre- 
vious example. Specifically we apply spherical boundary 
conditions (sbc), where the dynamic filaments are con- 
fined to a sphere of radius R. For comparison we also 
consider periodic boundary conditions (pbc). 
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FIG. 4. Radial component of the filament field F(r) for 
default parameters in 3D. For pbc there is no bias in the 
network structure as expected due to translational symme- 
try. For confining sbc, one observes approximately a linear 
behavior as predicted by the analytical theory. 
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FIG. 5. Distributions of plus- (black) and minus-particles 
(orange) for spherical boundary conditions. (a) Snapshot of 
particle configurations (Slice plane), (b) Radial density dis- 
tribution of vesicles. 



Due to translational symmetry for pbc and non-biased 
dynamics there is supposed to be no overall bias and the 
field F vanishes. This is reproduced by simulations as can 
be seen in Fig. 4. However, if the filaments are confined 
to a sphere, the theory predicts that F has the same form 
as an electrostatic field inside a homogeneously charged 
sphere. For that geometry, Gauss' divergence theorem 
yields 



F(x) 



for r < R . 



(4) 



This linear behavior of F, with average radial orienta- 
tion of filament plus ends towards boundaries, is indeed 
reproduced by simulations of filament ensembles (Fig. 4). 

Particles follow the emerging radial bias of the fila- 
ments leading to a separation of particle species (Fig. 5). 
Therefore, in the case of spherical boundary conditions, 
plus-particles accumulate at the boundaries while minus- 
particles are depleted at the boundaries and have a slight 
tendency towards the center. Analogue to electrostatics 
we conclude that, for an arbitrary geometry, random nu- 
cleation dynamics leads, in general, to an average bias of 
plus ends towards the boundaries. 

However, for MT dynamics the Green's function needs 
some modifications. Since there is a chance for plus-end 
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FIG. 6. Illustration of the re-organization of the microtubule 
network after transection. After dissolution of filaments, they 
re-assemble randomly in the right part (yellow region). Ac- 
cording to our theory, plus ends of MTs (black arrows) are 
pointed towards the boundaries of this region. Adapted from 
[3] 



depolymerization, while minus ends are fixed, the Green's 
function has a finite range scale, given by the mean fila- 
ment length. For larger distances, the effect of a nucleat- 
ing filament is 'screened'. This means that the bias effect 
is only strong within a distance from the boundary that 
corresponds to the average filament length. Nonethe- 
less, the direction of the bias remains unaffected by the 
screening. For our model geometry, the mean filament 
length (80 r p ) is larger than the cylinder radius (40 r p ), 
such that screening effects can be neglected. 

Within this theory, we can also explain the observed 
distribution of particles in case of a transected axon: Af- 
ter dissolution of MTs, random nucleation and growth 
induces a bias of plus ends outwards of the nucleating 
region, arranging (schematically) as in Fig. 6 (only right 
part (yellow) is nucleating region). Together with the 
surviving microtubules at the left end, pointing to the 
right, this yields a trap for plus particles. In addition, a 
trap at the right tip emerges due to MTs pointing to the 
right. Our model, however, underestimates the strength 
of the particle traps compared to the experiments (Figs. 
2 and 3). This discrepancy supports the interpretation 
given in [3] that interactions between filaments (which 
we have neglected) indeed increase orientational correla- 
tions to enhance structural inhomogeneities in the fila- 



ment networks [15] that stabilize and amplify the par- 
ticle trapping. The bias predicted by our model thus 
can trigger an interact ion- driven self-organization of the 
filaments towards the observed patterns. 

To summarize, we studied a model for active trans- 
port on filament networks that evolve by isotropic and 
spatially homogeneous stochastic filament dynamics. We 
showed that if filament dynamics is confined within a fi- 
nite volume, the filaments become oriented following the 
geometry of the system, despite the absence of filament- 
filament interactions and chemical gradients. Applied to 
the geometry of a transected axon, this random dynam- 
ics is sufficient to explain the structure of experimentally 
observed vesicle traps, which is essential for axonal re- 
generation. 

We developed a linear theory which is analogous to 
electrostatics and describes correctly the alignment of fil- 
ament orientations. This theory is assumed to be generic 
for arbitrary geometries and predicts a preferred orien- 
tation of filament plus-ends towards boundaries. Hence 
we have shown that in confined geometries, concentra- 
tion gradients can generally be induced by randomly 
generated networks of polarized filaments, even without 
interaction-driven self-organization or external gradients. 

Our results indicate that the geometry of the cell may 
have strong influence on the development of the struc- 
ture of the filament network: Although our model re- 
sults do not reproduce the strength of vesicle traps, 
their structure is well described. Therefore we conclude 
that the mechanism presented in this work may select 
the emerging structures of interacting filament systems 
[7, 9, 10, 16]. This assumption is also supported by ex- 
periments providing evidence that the cell shape may also 
control the organization of the cytoskeleton [17]. 
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Marc Neef for fruitful discussions and Maren Westkott for 
image analysis. We also thank DFG Grants GK1276/1 
and GZ SA864/3-1 for financial support. P.G. is post- 
doctoral fellow of DA AD. 
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Supplementary Material 



I. FILAMENT AND PARTICLE DYNAMICS 

Tables I and II are summarizing the elementary steps 
of the particle and filament dynamics. Note that we dis- 
tinguish between microtubule (MT) and actin dynamics 
in our model in the following way [1]: In contrast to MTs, 
actin filaments can uncap and depolymerize at the minus 
end to perform treadmilling. MTs do not treadmill but, 
on the other hand, can perform dynamic instability at 
the plus end. They can turn to a plus-end depolymeriz- 
ing state, which can be rescued to turn to a polymerizing 
state. 

For all simulations we used a stochastic random se- 
quential update scheme and continuous space coordi- 
nates. 

II. CHOICE OF DEFAULT PARAMETERS 

The default parameters of our model were chosen in 
accordance with available experimental results or, if not 
directly accessible, according to established models. In 
Table III the parameters for microtubule and vesicle dy- 
namics are given. For dynamics of actin dynamics, we 
refer to the work [2], Table 4. 

The measured quantities in Table III (column 4) are 
displayed according to their sources referenced in column 
3. These measured quantities are transferred to the 
time and length scales used in the simulations (column 5). 

Comments: 

(i) The maximal tubulin density is reached, if all mi- 
crotubules are dissolved. Assuming 26 tubulin 
dimers per d s = 16 nm microtubule segments [1] 
(16nm are two twists of the microtubule structure), 
20 microtubules per jim 2 [3] and free tubulin con- 
centration of in healthy axons [4], one obtains 
59/iM. 

(ii) All the minus ends of microtubules are bound to 
7— tubulin which also serves as nucleation seed. 
Assuming an average microtubule length of 4 fim 
[5] and 20 microtubules//im 2 [3] one obtains 5 
7— tubulin nucleation seeds per /im 3 . 

(iii) All nucleation seeds are already present in the be- 
ginning. Each microtubule is hence initiated by a 
growth process of a bare 7— tubulin. 

(iv) In reference [6] the binding rate for kinesin in de- 
pendence of the distance x of the binding filament 
is given by 3^exp(— x 2 /2cr 2 ) with n = 175s -1 and 
o~ = 6nm. The average binding rate within a dis- 
tance alb — 80nm from a filament of length L hence 

is uj a = l/(7rd 2 b L) J L J Q db 27rre"^ dr dz = 5.9s' 1 . 



The integration space is due to the cylindrical form 
of the binding region. 

(v) We adjusted the entry and exit rates of particles 
such that the observed density in the healthy axon 
was approximately as in the referenced work. 

(vi) The length of the microtubule depletion zone in 
the experiments [7] was about 100 jim. Due to 
computational constraints we were able to simulate 
a system of length L = 200 lu=2Q /im, giving the 
length of the depletion zone L/2 = 100 lu. Adjust- 
ing to the proportions of the experimental picture, 
we chose a cylinder radius 20 lu. However, we also 
tested other system sizes which did not affect the 
generic structure of the particle traps. 



III. TIME EVOLUTION OF THE 
SIMULATIONS: TRANSECTED AXON 

Here we show some intermediate time configurations of 
particles and filaments both in the case of a plain cylin- 
drical geometry (with a round tip) and a cylinder with a 
tapered region, mimicking the geometry of the axon more 
faithfully. One observes that in both cases accumulations 
of plus particles occur left of the depletion area and at 
the tip, similar to the experiments. However for the plain 
cylindrical geometry it takes longer to achieve this state. 

One also observes that in the beginning, the parallel 
bundles are also present in the right part of the system, 
while accumulation of vesicles already begins. Hence the 
accumulation is not due to a diffusive bottleneck in the 
middle of the system. 
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Process Particle state (s) 



Description 



Parameter name 



Diffusion D 
Step A 
Attachment D— )>A 



Detachment 



A^D 



Detached particles move in a random direction. Step 
widths are uniformly distributed between and 21 d ■ 

Attached particles move to adjacent subunit in 
(+) -direction. 

Particles bind to subunits if their distance is less than 
db, becoming 'attached'. 

Particles detach from filament. 



TABLE I. Particle dynamics. A='attached'; D='detached'. 



Process 



Description 



Probability 



Nucleation New filament created; arbitrary direction. The process 
needs nucleations seeds (e.g. 7-tubulin for microtubules) 
and monomers. (— )-Cap inhibits shrinking. 

Growth New sub-segment added at (+)-end. 

Shrinking Segment removed at ( — )-end if not ( — )-capped or at 

(+)-end if in shrink state. 
Uncapping ( — )-Cap is removed allowing depolymerization hence 

(only actin). 

Catastrophe (+)-end turns into shrink state, (only MT) 

Rescue (+)-end turns into growth state again, (only MT). 



Ldn Pmon Pnuc 

pnuc — nucl. seeds 

pmon = density of free monomers 

Pmon 

cop- 



Meat 
0J res 



TABLE II. Filament dynamics. 



FIG. 1. Configurations of particles and microtubules for simulations of the model at distinct times, (a) Geometry of a transected 
axon, (a.l) t=200 (a.2) t=3000 (a.3) t=5000 (a.4) t=24000. (b) Cylindrical geometry, (b.l) t=200 (b.2) t=5000 (b.3) t=20000 
(b.4) t =32000. The qualitative structure does not depend on the details of the geometry. 
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Parameter name 


Reference 


Reference Value 




Model parameters 


Filament dynamics: 


growth rate uj 3 


[4] 


1.5 \1m\1M 1 min 1 




0.015 lu 3 tu- 1 


shrink rate uj s 


[4] 


50pm min~ x 




O.Stu- 1 


microtubule density pMT 




20 pmT 2 




0.2 lu~ 2 


max. monomer density pmon 


7-tubulin density [1, 3, 4] 


59 pM ( See comment (i)) 




35lu~ 3 


7— tubulin density p 1 


[3, 5] 


5 pm~ 3 (See comment (ii)) 




0.005 lu~ 3 


nucleation rate oj n 


See comment (iii) 


w g p 1 ptub 




see above 


catastrophe rate oj c 


[8] 


0.014s- 1 




0.0001 tii' 1 


rescue rate oj r 


[8] 


0.04 s- 1 




0.0004 tii' 1 


Particle dynamics: 


particle radius r p 


[9] 


50nm (average) 




0.5 lu 


binding distance db 


[1] 


80 nm (length of kinesin) 




0.8 lu 


subunit distance d s 


[1] 


16nm (1 stepping period of kinesin) 




0.16 lu 


step rate p 


[10] 


velocity — 1.3 pms -1 




O.Utu' 1 


attachment uo a 


[6] 


5.9 s -1 (See comment (iv)) 




0.059 tu- 1 


detachment ujd 


[10] 


run length « 7 pm 




0.0018 tu- 1 


diffusive step length Id 


[11] 


diff. const. D = 2.5 x 10 -10 cm 2 s' 1 




\/2Dtu = VoMlu 


particle density p® 


[12] 


7.6 pm~ 3 




0.0076 lu~ 3 


entry rates a 1,2 


[12] 


See comment (v) 


ai 


= 0.15^-\a 2 = Q.012^- 1 


exit rates f3\,2 


[12] 


See comment (v) 


Pi 


= 0.003 tu~ 1 ,p 2 = 0.2 tu- 1 


System size 


See comment (vi) 


See comment (vi) 


length L = 200 lu, radius R = 20 lu 


MT depletion zone 


See comment (vi) 


See comment (vi) 




L/2 



TABLE III. Default parameters of the model which are biologically motivated by transport of vesicles in axons. The referenced 
values are either based on experimental data or existing models for microtubule dynamics [4]. Model parameters are chosen 
to be in the order of magnitude of referenced values, fitted to time and space scale of the simulations. Length scale: 1 lu = 
lOOnm = 2r p lpM = 0.6lu~ 3 . Time scale: 1 tu = At = 0.01s which is one time step in the simulations. We consider a 
cylindrically shaped system length L = 200 lu and a radius of L/10. In the transected state, it has a spherically shaped tip and 
a neck region right of the middle with a minimum radius 0.6L. 



